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Abstract 

The pseudogap is a key property of the cuprate superconductors, whose understanding should 
illuminate the pairing mechanism. Recent experimental data support a close connection between 
the pseudogap and an oxygen-driven C4 symmetry breaking within the Cu02 plane unit cell. Using 
ab initio Molecular Dynamics, we demonstrate the existence of a strong nonlinear electron-oxygen 
vibrator coupling in two cuprates. In a mean field approach applied to this coupling within a model 
Hamiltonian, we derive a C4 splitting/pseudogap phase diagram in agreement with experiment - 
providing an explanation for the pseudogap phenomenon from first principles. The implications 
for superconductivity and the Fermi surface arc effect are discussed. 



In 20 years no consensus has emerged on the mechanism for <i-wave high temperature 
superconductivity (HTS) in cuprates or for their characteristic pseudogap (PG) Al- 
ternatives such as the electron-lattice coupling which mediates pairing in conventional su- 

n n n 

perconductors [2j, and purely electronic mechanisms [3,14], have been widely discussed. A 
possibly decisive input comes from recent work showing a strong correlation between 
the occurrence of the PG in the underdoped region and the breaking of C4 symmetry (i.e. 
the a and b directions become nonequivalent) in the planar oxygen vibrational amplitudes. 
This observation points to the origin of the PG as a lattice instability, and is a very im- 
portant clue - consistent with evidence such as the doping-dependent oxyg en isotope shift 
7, s| and the HTS-induced softening in oxygen vibration frequency 9), [lo| - that electron- 
lattice coupling is the mechanism for HTS. Paradoxically, ab initio calculations do not find 
a strong conventional linear electron-lattice coupling in the cuprates [ll|, 12 ] , leading to the 
proposition that a primarily nonlinear coupling lies at the heart of HTS. 

In this paper we use ab initio molecular dynamics (AIMD) to show that electron coupling 
to the vibrations of O in the Cu02 plane is in fact very strong - but it is second order in the 
oxygen displacement, in contradiction with the assumptions underlying the earlier electron- 



lattice coupling work 
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121 ] . To explain why the coupling should be second order, we argue 



from symmetry that the local effect on the electrons of a radial displacement of oxygen in the 
Cu-O-Cu bond should be independent of the displacement's sign, and hence second order 
in the O displacement. Based on this symmetry argument, two of us earlier put forward 
the Fluctuating Bond Model of HTS [131] which is empirically based on the idea of second 
order coupling, an idea which had been proposed earlier {14I 



16j. The FBM successfully 



predicted C4 symmetry breaking 13j and enabled superconducting properties such as d- 
wave symmetry [17 1 to be understood. Unlike conventional superconductivity, the FBM has 



not yet been parametrized from first principles, and here we show that t 



FBM coupling energy does indeed have the value required in the FBM 
key assumption of a giant nonlinear electron-lattice coupling. 



re value of the 



131 ]. validating its 



The technique of ab initio molecular dynamics (AIMD) 18] solves the ionic equations of 



motion on a Born-Oppenheimer potential energy surface obtained by solution of the many- 
electron Schrodinger equation in local density approximation, a technique well suited to 
the present problem in that it avoids the constraint of a linearized interaction. We start 
by applying AIMD to the C4 symmetry-breaking in two cuprates. First we consider the 
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FIG. 1: ab initio Molecular Dynamics calculation at T = 4K of the structure of metallic La2Cu04 
(214) [19(, blue spheres, La, green spheres, Cu, red spheres, O. a, Undistorted setup structure, b, 
Equilibrated structure showing vertical displacements of the planar oxygens correspond to rotations 
of CuO"6 octahedra about alternate x- and y- axes in planes stacked along the c-axis - the LTT 
structure found at low temperature in metallic 214 phases. 



well-established Low Temperature Tetragonal (LTT) structure [19( found in underdoped 
metallic 214 materials, which was early on linked to non-linear electron lattice coupling 
[20! . Figure 1 shows the LTT structure is predicted correctly by AIMD for the metallic 
La2CuC>4 at T = 4K. All Cu02 planes undergo C4 symmetry-breaking in which one half 
of the oxygens undergo a z-displacement, the n/2 rotation of successive planes around the 
c-axis ensures overall tetragonal symmetry. The symmetry breaking by ^-displacement in 
the 214 family contrasts with the symmetry breaking by xy-displacement in the 2212 and 
oxychloride (Ca 2 Cu0 2 Cl 2 ) materials 6]. 

We now establish a model-independent global connection between oxygen PE surface and 
the number of electrons by extending the Fig. 1 type of calculation to a range of dopings. 
The oxychloride system Ca2- a Naa;Cu02Cl2 is computationally attractive and doping can be 

culated oxygen PE surfaces being shown in Fig. 



controlled via the Na fraction x, the c&. 

2. AIMD accurately reflects experiment (| in that the oxygen displacemements are in the 
xy-plane (Fig. 2 inset) (but cannot capture its non-Born-Oppenheimer features). The PE 
surfaces in Fig. 2 are strongly dependent on doping, their minimum is off the bond axis 
(leading to the Fig. 2 inset distortion) at low doping, transitioning to an on-axis minimum 
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FIG. 2: PE curves V(u) for oxy chloride material as a function of oxygen in-plane distortion u 
(see insert) for different dopings (see labelling on curves). Doping is implemented by fractional 
substitution of Na for Ca. 

(stability of oxygen on the bond axis) at high doping, and can be parametrized in the form 



Here u is oxygen displacement from the bond axis, the zero-doping bond force constant x is 
negative (oxygen unstable at zero doping), and the electron-lattice coupling V is positive, 
representing stabilization of the unstable Cu-O-Cu bond with increasing doping p. The 
positive quartic term w represents the confinement of oxygen in the local lattice cage. The 
results of Fig. 2 and Eq. (pQ) clearly show that the oxygen potential energy surface is strongly 
dependent on global electron number, with a coupling term in V going as the square of the 
oxygen displacement, as was argued above. If the Cu-O-Cu bond is sensitive to the global 
number of electrons, we can also anticipate a strong response of the oxygen to the local 
electron number. 

To calculate the local form of electron-lattice coupling requires a model. We start from 
the Emery tight-binding model based on Cu ?>d x i_ y i and the oxygen 2p x /2p y orbitals that 
have a symmetry in the Cu-O-Cu bond. The coupling is determined by considering how 
the matrix elements vary with oxygen displacement, and projecting the result back into a 




(1) 
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one-band basis of just the 3d x 2_ y 2 orbitals (see Supp. Mat.). The resulting complete FBM 
Hamiltonian Eq. ([2]) is written in mixed representation electronically. It considers only one 
vibrational mode per oxygen for simplicity: 



k,cr (ij) 



iJ U ij- 



Equation (J2J) embodies the conventional 1-band model [2l| with a single 2>d x i_ y i orbital per 
Cu atom located at site % in the Cu0 2 plane. Here cf a (cj >(7 ) is the creation (destruction) 
operator for the 3d x 2_ y 2 orbital on site i of spin a, with number operators n i)(J = cf a Ci t(7 , 
and c£ CT (ck, CT ) are the corresponding operators for the band states of wavevector k, etc.. 
A sum over (ij) implies that each nearest-neighbor bond ij appears only once in the sum. 
The variables Pij, Uij are the conjugate momentum and position coordinates of oxygen 
in bond ij, M is oxygen mass and Xo the bare oxygen force constant, w is the quartic 
interaction, and v a (positive) electron- vibrator coupling constant. The operator iV^ = 
J2a ( n i°- + n j° ~ c t,a c j,o- ~ c ^a c i,o) represents the number of electrons in the r/'-antibonding 
orbital [cj )(7 — Cj >a ] / y2. 

In more detail the first term of Eq.(j2]) is the electronic band energy, and the second 
term represents the harmonic part of the oxygen Hamiltonian, to which is added the third 
term, a quartic interaction which confines the oxygen and is needed to fit the PE curves 
in Fig 2. The nonlinear electron- vibrator coupling is represented in the fourth term, which 
involves the vibrator coordinate quadratically, and the electronic operator Nij. Dumping 
electrons iVy into the antibonding orbital in bond ij softens the vibrator ij, as implied by 
Eq. (CQ) and by the trend in Fig. 2 - it is also reasonable on intuitive chemical grounds. Note 
that the c^ a Cj ia + Cj ff Ci,o term in N^, which originates from processes where an electron hops 



from Cu i to O and then to Cu j, is identical to the original FBM [131 ]. while the present 
analysis adds a new term — {rii a + originating from processes where an electron hops 
from Cu to O and then ends up back on the same Cu. 

The FBM coupling is strong between the oxygens and the saddle points at X= (tt, 0) and 
Y= (0, 7r) in the band structure - the energies at X and Y are normally degenerate, but the 
degeneracy is split by C4 symmetry-breaking. The bare electron-lattice coupling constant 
v in Eq. (T5J) can be simply determined by displacing the x-oxygens and determining the 



shift in the band structure eigenvalue at X. Any effect of a chemical potential shift can be 
removed by displacing the ^/-oxygens and subtracting the two effects (see Supp. Mat.). The 
results are collected in Table I for the 214 and oxychloride materials. The values of the 
quartic interaction w are obtained from the quartic coefficient of the fit to the Fig. 2 curves, 
and similar ones for the 214 material. It is found that the coupling v is relatively small 
for vibrational polarization along the Cu-O-Cu bond, so we only considered polarizations 
transverse to the bond in Table I. 



Table I: FBM Parameters for 214 and Oxychloride Materials 



Polarization 


v 2 ia (au) 


v oxy (au) 


w 2 u (au) 


^oxy (au) 


xy _L to bond 


0.0163 


0.0182 


0.053 


0.090 


z J_ to bond 


0.0174 


0.0202 


0.122 


0.106 



Note that the coupling strengths v in Table I are rather material- and polarization- 
independent, while the anharmonic coefficents w vary by less than 2.5 x. 

The pairing energy in the FBM is second order in the coupling v, with the functional 
form [13J v 2 /w. In the current form Eq.(j2J) it is K' = Av 2 /w (see Supp. Mat.). From 
Table I K' can be estimated to be about K' ~ 0.5 eV, matching the giant magnitude 
required [3] for high temperature superconductivity. This is the first meaningful ab initio 
prediction of the electron-lattice coupling constant in HTS, and points the way to an eventual 
quantitative theory of HTS by extending the existing body of work describing conventional 
superconductors. 

Here, as an application of the model Eq.(j2]), rather than again presenting results for the 
gap equation [l3j which will be done elsewhere, we present an important new result, the 
pseudogap phase diagram. This is done by the mean field decoupling of Eq.(j2J) (see Supp. 
Mat.), (a) NijU 2 j — > (iV^-uf- , leading to a softening of the vibrator proportional to the 
number of antibonding electrons (Nj), and (b) Nju^ — > N^ (u 2 ,) wherein the nearest- 
neighbor hopping (c^ a Cj i(T + dj Ci >(T ) is reduced by increased vibrator amplitude (u 2 ^. In 
this mean field approximation the model decouples into exactly soluble band structure and 
anharmonic oscillator problems, linked by the mean field quantities (N^) and (ufj), which 
are solved for self-consistently assuming they are uniform in space. 

When the expectation values (Nij) and (u 2 j) in the x- and ^/-directed bonds are found to 
differ, C4 symmetry breaking in both vibrator amplitudes and the electronic structure occurs. 
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Symmetry breaking in the oxygen vibrator amplitudes is exactly the effect measured in the 
STM experiments of Ref. {f| (consistently we predict that the higher/lower-amplitude oxy- 
gens in the Ref. |f| R-plots should be dark/light streaks (see Supp. Mat.), as observed). C4 
symmetry breaking in the electronic structure, wherein the degeneracy of the saddle points 
at X= (tt, 0) and Y= (0, tt) is split, leads to a <i-type PG A ps (k) ~ A ps (cosk x — cosk y ) /2 
(see Supp. Mat.), where the PG, A ps , can be positive or negative in sign. The PG A ps is 
one half the splitting between the energies at X and Y. The work of Ref. [5] shows that the 
experimental PG is intimately associated with the C4 splitting - confirming the correctness 
of identifying A ps with the experimental pseudogap. The predicted PG phase diagram is 
shown in Fig. 3, and is seen to reproduce the known experimental phase diagram for the 



pseudogap 



24 



25] . At low temperature the pseudogap ranges from a maximum of 



about lOOmeV on the underdoped side, decreasing with increasing doping |5J, |2i 
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FIG. 3: Contour map of pseudogap in temperature/doping plane showing decrease with doping, 
and with temperature, until it vanishes at phase boundary T*. Contours labelled by pseudogap 
A ps in intervals of 13.75 meV. For experimental A ps magnitudes see Ref. j^ . 

Note that experimentally, the C4 splitting is characterized by a nanoscale domain struc- 
ture 6|. In our picture this implies that the sign of A ps will vary spatially, in a manner such 
as A ps (k) ~ cos^Aps (cosfc^ — cosk y ) /2, where 6 is some local phase determined by the 



local nanoscopic disorder. The Fermi surface in the C4 split phase will be sensitive to the 
sign of A ps . The inset in Fig. 4 shows the two FS corresponding to the extreme cases 8 = 
and 9 = it, and it is reasonable to assume that a nanoscopically varying order parameter 
will lead to FS smearing as illustrated by the shaded region in the Fig. 4 inset. 

The smearing of the FS by the PG, seen in the Fig. 4 inset, is seen to be least near 
the nodal line k x = k y , a manifestation of the <i-type nature of the PG A ps (k), and to 
increase as one goes towards the SP's at X and Y. If we make a measurement on some 
energy scale E, it might be expected that the FS will be well-defined at k-points where the 
local PG is less than E, A ps (k) < E, but that it will appear smeared out on energy scales 
where the local PG is larger than E, A ps (k) > E. As an example of this phenomenon we 
could take E to be temperature, E = k^T . According to this argument, there should then 
be a boundary between the resolvable and unresolvable sections of the FS arc defined by 
A ps (k) = k B T. In fact this FS arc effect [27| has already been observed a couple of years 
ago, as shown by a comparison of the data with our heuristic model in Fig. 4, which is 
in good agreement with experiment 28] (this approach does not include the temperature- 
dependent quasiparticle lifetime broadening {2^ , which we hope to derive in a future paper 
dealing with quasiparticle dynamics). 

The approach here has been to focus on the electron-lattice interaction while ignoring 
the strong electron correlation effects. It would be productive to consider the nonlinear 
electron-lattice coupling in the context of a strongly correlated (large-?/) system, where the 
modulation of the Cu-Cu bond by oxygen vibration is expected to be as strong as in the 
noninteracting case, and the phase diagram can be extended to include magnetic effects. 

In conclusion, using AIMD we have demonstrated the instability of the Cu02 plane 
oxygens in three cuprates which relates closely to experiments such as the LTT structure 
in the 214 family and the C4 symmetry breaking in the 2212 and oxychloride materials. 
The nonlinear electron-lattice coupling deduced from these calculations is large with the 
magnitude required by the FBM theory. Using our deduced Hamiltonian Eq. (2), we 
are able to derive the pseudogap phase diagram, harmonizing the pseudogap data with the 
observed C4 splitting, and to understand the temperature-dependent Fermi arc observations. 
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I. AB INITIO MOLECULAR DYNAMICS 

The Car-Parrinello ab initio MD [1] studies presented herein were performed using Kohn- 
Sham Density Functional Theory. The CPMD approach is based on the Born-Oppenheimer 
approximation. It describes the metallic state, but not the magnetic, insulating state, of the 
cuprate materials. The computations were checked by comparison with the band structures 
of copper and the 214 material [2|. The Fig. 1 results were obtained by starting with the 
Fig. la configuration and running to equilibrium. The Fig. 2b results were obtained by 
constraining the oxygens in the manner shown in Fig. la at a specified value of the oxygen 
displacement, and running to equilibrium - they therefore represent " ion- relaxed" potential 
energy curves. Doping was achieved by varying the K/Ca ratio in the ensemble. Samples of 
4x4 and 6x6 unit cells were run in periodic boundary conditions. 



II. PROJECTION OF 3-BAND MODEL ONTO 1-BAND MODEL 



In matrix notation consider a <i-subspace and a p-subspace, represented by the Hamilto- 
nians H d and H p respectively, connected by the coupling matrix V pd , the Hamiltonian then 
being 



H 



H d V dp 

ypd H p 

Projecting onto the <i-subspace in perturbation theory 



(1) 



H d = H d + V dp (e d - H P Y l V pd . 
if i, j are <i-sites, and I, m are p-orbitals 

H*v = eAi + £ V* (e d - H*)£ V*. 

Lm 



(2) 



(3) 



Now we shall neglect the pp hopping matrix elements (Emery model), when I = m and 
the V^'s are nearest-neighbor hopping matrix elements defined as t p d > (see Fig. 1). There 
are 2 processes, 



1. i,j nearest neighbor on the d-lattice, when the 2 V's have opposite sign (Fig. 1) 



2 






2p x 



FIG. 1: 2p x , 2py and 3d x 2_ y 2 orbitals in Cu0"2 plane, illustrating 2p to 3d hopping integral t p d- 

2. % = j, when the 2 V"'s have same sign 
giving 

H d = j2 wt* fa* + n i°) - E — x *v ( 4 ) 

where e p d = — e p > is the "oxide gap" between the oxygen 2p orbital energy and the 
higher- lying Cu 3^2.^2 orbital energy, a is spin, p^ is the p-orbital between ci-sites i and j, 
and the bond order operator is 

a 

Let us assume that the oxygen motion in some direction is x, and that it enters the 
3-band hamiltonian via the pd hopping integral 

tpd -> t pd - v pd x 2 , where v pd > 0, (6) 
then to order v p d, and defining t = t^ d /e p d 

= (e d + 2t) E n - - t E X V ( ? ) 
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Restoring our original notation [3|] 2t p dV p d/e p d = v/2y/nn s , when the coupling v is seen 
to be positive 

H d = {e d + 2t) nia - t ^ X l3 (8) 

(hi) 

V \ -v , , o f \ -v „ 2 



5^ K„ + n ja ) x% + — = Yl X i 



X: 



2Jnn s ' J 2Jnn s ' J 

We retrieve our previous 1-band model (next-nearest and next-next-nearest neighbor 
hoppings are dropped due to neglect of £ pp ), but with an extra term diagonal in rf-space. 
As regards the vibrator, the effect of the new term is to weaken the oxygen parabolic 
potential linearly with increasing electron occupation of the band (or stiffen the vibrator 
with increasing hole occupation). The number operator term is dominant over the hopping 
term (<X> maximizes at ~ 0.6). 

Let us now alternatively assume that the oxygen motion enters the 3-band hamiltonian 
through the interaction of the electrostatic potential with the charge on the oxygen 

e pd ^> e pd + v p x 2 ; (9) 

where v p depends on a Madelung sum. In an ionic crystal it is arguable that the sign of v p 
will be positive since the environment of a negative ion typically consists of positive ions, 
so as the O-ion approaches them the local oxide gap e p d becomes larger. However in a 
perovskite structure the issue needs specific calculation. 
Expanding to first order 

I I v p x 2 (1Q) 



£pd V p X d p d £ p d 

Returning to Eq. (j4j), we insert the foregoing expansion into the 2 terms to obtain 

AH* - £ (n ia + n jc ) x% + ^ £ X ijX %. (11) 

The effect of the oscillator correction ffTTj) from this mechanism can be absorbed into (jSJ), 
giving the same final result jS]) but with 

v 

9 , = (2t pd Vpd + tv p ) I e pd . (12) 

The sign of v will be positive if the t p dV p d term in parenthesis is dominant, or if v p is positive 
as argued above. 

In this section we have formally derived the FBM coupling, showing the approximations 
involved explicitly, and demonstrated the existence of a new term in the coupling. 



III. FBM HAMILTONIAN 



The FBM Hamiltonian involves three pieces 

H = H v + H e + H ev . (13) 

In H the Cu sites, which define the unit cell, are defined as 2D integral-component vectors 
i = (i x ,i y ) (lattice constant is taken as unity). The two oxygens in each unit cell i are 
located at the sites i+a/2, where a is a unit vector along the x- or y- axes, hence a defines 
whether the oxygen is in a Cu-O-Cu bond oriented along the x- or y- direction. 

In the vibrator piece H v the oxygen degree of freedom is an n-component vector Xi + s/ 2 , 
where n — 1 if a single mode is dominant (as assumed in the manuscript), n = 2 if the two 
modes transverse to the Cu-O-Cu bond are roughly equivalent, or in a case now considered 
unlikely (as the along-bond mode is found to be weakly coupled) n = 3 if the two transverse 
modes and the along-bond mode can all be considered equivalent. H v is given by 

r 1 2 , XO 2 W 



(14) 



In H v the scalar products Xi + a/ 2 ■ Xi+a/2 are abbreviated to x? +s ^ 2 , and a momentum Pi+s/2 
conjugate to coordinate x i+ s/ 2 is introduced, to define the vibrator kinetic energy, with m 
the oxygen mass (M in the Ms.). The "bare" bond force constant is xo- The quartic term, 
with coefficient w, is assumed in the degenerate case to be radially (n = 2) or spherically 
(n = 3) symmetric. 

The electronic piece H e is 

H e = -lY^t(i-])clc^ (15) 

where cf a (ci y0 ) denote respectively the creation (destruction) operators for the ?>d x 2_ y 2 orbital 
(or, more rigorously, the d^.^-type Cu3<i-02p antibonding Wannier function) on lattice site 
i of spin a. The strongest interaction is the nearest neighbor hopping integral £(±1,0) = 
t(0, ±1) = t, (t is positive), followed by the next-nearest neighbor interaction t(±l, ±1) = t', 
(t' is negative) and then the 3rd-nearest neighbor interaction t(±2, 0) = t(0, ±2) = t" (t" is 
positive). The band eigenvalues of (I15p are 

ek = — 2t(cos k x + cos k y ) — 4i' cos k x cos k y — 2t"(cos2k x + cos 2k y ). (16) 



The model band structure has a minimum at T (k = (0, 0)), a maximum at Z (k = (jr, 7r)), 
and saddle points (SP) at X (k = (0,7r)), and Y (k = (tt,0)). As a result of the saddle 
points, located at e§p = 4t' — 4t", the density of states (DOS) has a logarithmic peak (van 
Hove singularity or vHs) at esp which is found from ARPES and band structure calculations 



the resulting high 



for near-optimally doped systems to lie close to the Fermi level 
DOS at the Fermi level strongly enhances the FBM coupling. The total band width is 8t. 
The electron- vibrator coupling piece is 

y 



H e 



a/2 



X\ + a/2 — ^2 ( C tv Ci+S ' CT + C t+S,a C h° 



(17) 
(18) 



where the bond order operator X is associated with the oxygen site at the bond center, 
and we have defined in the mixed degeneracy factor {nn s ) 1 , where n s = 2 is the spin 
degeneracy, to make the term of order y/nn s , motivated by a version of large-iV theory 
jointly expanding in 1/n and l/n s . In Ref. 3[ only the X-piece of (JTTj) was included. 

The combination — [n- ha + n. &>a ) + Xi + s/2 can also be written in more compact form, 
defining the antibonding orbital \a, i+a/2) = (|i) — |i+a)) /\/2, with number operator n? +s ^ 2 
(summing over spin). In the manuscript n a is simplified to N . 



n 



i+a,aj 



Xi+5/2 



(19) 



The complete Hamiltonian H = H v + H e + H ev is then 
y 



H 



E 



1 2 XO 2 W ( 2 \2 



1 - JJ Ci )CT Cj, CT 



(20) 



+ 



'nn q 



+S/2 n i+S/2- 



i,o=x 



Note that in Eq. (l2"U|) if = w 2 /u; defines a coupling energy. 



IV. DETERMINATION OF COUPLING v 



Calculation of the oxygen PE surface as a function of doping is not an ideal approach to 
calculating the FBM coupling constant for two reasons. The PE surface calculations such 
as illustrated in Fig. 2 in the Ms. are relaxed surfaces, i.e. when certain oxygen coordinates 
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are held fixed all other atoms are allowed to find their minumum energy. This differs from 
the "vertical" couplings entering the FBM Hamiltonian, where atoms other than planar 
oxygen are assumed fixed in place. Secondly the coupling in the FBM Hamiltonian is to the 
number of electrons n a in the antibonding orbital, which mainly involves states at the top 
of the rf-band and will be filled mainly by adding electrons rather than holes as was done 
(for reasons of computational stability) in Fig. 2 of the Ms.. 

The method adopted to calculate the coupling strength v is based on comparing the 
shift in band structure energies when the oxygen location is perturbed with the same shift 
deduced from the FBM Hamiltonian. The FBM coupling (third term in Eq. (1201) ) leads to 
splittings in the tight-binding band structure. If all oxygens in the x-oriented bonds are 
globally shifted by u x , and all oxygens in the ^/-oriented bonds by u y , there is a splitting of 
the band energy between the band energy ex at the saddle point (SP) X= (7r,0), and ey 
at Y= (0,7r), given by ex — ey = \/2/nv (u 2 x — u^j . By numerically calculating the band 
structure with first the x-oxygens displaced, and then the y-oxygens, and subtracting the 
corresponding band structure energies energies at, say, the SP X, any isotropic shift resulting 
from displacing a single oxygen can be cancelled out and the coupling v determined. The 
results are shown in Table I. 

V. MEAN FIELD APPROXIMATION 

Mean field theory is a useful step in investigating the properties of many models. In the 
FBM, the mean field approximation decouples the electronic and vibrational parts of the 
Hamiltonian. In the vibrational part, an expectation value of the electronic terms shifts the 
oscillator harmonic frequency, the expectation value being assumed spatially uniform, but 
it can be different in the x- and y- bonds (in this section we return to the notation in the 
Ms.): 




(21) 



v 



{2-2p+{cf r c jiff + c+ a c i) , 




+ 



2^2 
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H mb can easily be diagonalized in a harmonic oscillator basis. In the electronic part, the 
expectation value of the square of the oscillator amplitude has been taken, 

H el = WK* + -£j= ]T l c t c ^ + c t^A «) , ( 22 ) 

k,cr V (ij) t a 

giving a band structure problem in which there are new nearest-neighbor hopping terms 
(v/2y/2) [cf Cj ta + dj Ci t(T \ (the uniform shift represented by the number operator 

terms does not change the band structure and is omitted) with the effect of reducing the 
nearest-neighbor hopping integral. Allowing the oscillator amplitude squared for the re- 
directed \Uij) x and y-directed (u^-) bonds to be unequal (the C4 symmetry-split case), the 
band structure is changed to 

? k = e k + -j= (u 2 ij ) x cos k* + -j= { u %) y cos k y (23) 

Using the band structure ?k ( 1231) the expectation values (cl Cj >a + Ci >er ) for x-oriented 
and ^/-oriented bonds are calculated, hence defining two quartic Hamiltonians pip whose ex- 
act solution yields the squared vibrator amplitudes (uij) x and • These interconnected 
electronic and quartic problems are then solved self-consistently as regards the expectation 
values. The parameters used were similar to Table I, v — 0.0198 au, w = 0.085 au, the 
oscillator bare force constant was xo — —0.0225 au. The band structure is parametrized by 
the (negative of the) hopping matrix elements, the nearest-neighbor hopping matrix element 
t = 0.25 eV, next-nearest-neighbor hopping m.e. t' = —0.05 eV, and third next-nearest- 
neighbor hopping m.e. t" = 27.2 meV. 

We can rewrite the effective band structure as 

?k = ?k + ^ A ps (cos k x - cos k y ) , (24) 

where A ps = (y{ u %) x '~ { u 1j) y ^ * s ^ ne pseudogap, and the renormalized nearest- 

neighbor hopping (d/2\/2) (y{ u fj) x + ( u ij)y) ^ s absorbed into e^. The experimental data 
{f| show that the pseudogap is not uniform over the sample as we have, for simplicity, 
assumed, but the coherence length over which the sign of A° varies is quite short, only a 
few lattice spacings. Probably as a result of this nanoscopic domain structure, the phase 
boundary of the pseudogap region is not typically found experimentally to constitute a true, 
sharp, phase boundary [3]. 

The variation of pseudogap with doping at low temperature seen in the contour plot 
(Manuscript Fig. 3) is similar to thet seen in experimental data (see Fig. 2) 
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Histograms of measured energy gaps A from a sequence of samples with 
different dopings, black being strongly overdoped and blue strongly underdoped. 
FIG. 2: Histograms of measured energy gaps A from a sequence of samples with different dopings, 

black being strongly overdoped and blue strongly underdoped 



VI. INTENSITY VARIATION IN EXPERIMENTAL i2-PLOTS 

In order to model the experimental behavior in the STM experiments [6] on C4 symmetry- 
split systems, we calculated the projected DOS for a 3-band model with the basis of oxygen 
2p x and 2p y orbitals and Cu 3d x 2_ y 2 orbitals shown in Fig. 1. The pda hopping matrix 
element is t p d = 1.12 eV. There are pp hopping matrix elements between nearest-neighbor 
2p x and 2p y orbitals given by t pp = —0.528 eV, and an oxide gap — e p = 6 eV. A spatially- 
uniform pseudogap is introduced by modifying the t p d matrix elements to t Pxd = t pd + At 
(i.e. for the lower vibrational amplitude oxygen) and t Py d = t pd — At (i.e. for the higher 
vibrational amplitude oxygen), where At = 0.0375 eV (the argument below only depends on 
these being semiquantitatively correct). The results for the DOS projected into the oxygen 
2p x orbitals (lying in x-oriented Cu-O-Cu bonds - see Fig. 1) and oxygen 2p y orbitals are 
different, as seen in Fig. 3. The DOS peak associated with the van Hove singularity is 
seen in Fig. 3 to be split, the peak above the Fermi level being localized only on the lower 
vibrational amplitude oxygen, and the peak below the Fermi level being being localized only 
the higher vibrational amplitude oxygen. The STM i?-map technique {f| for detecting the 
C4 splitting experimentally involves the ratio R of the tunneling current into the empty DOS 
to the hole current into the filled DOS. Evidently from Fig3 , R is predicted to be large on 
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FIG. 3: Oxygen projected 2p-DOS for oxygens in x-oriented and y-oriented bonds. The peak above 
the Fermi level is for the lower vibrational amplitude oxygen, and the peak below the Fermi level 
is for the higher vibrational amplitude oxygen. 



the low amplitude oxygens and small on the high amplitude oxygens, in agreement with the 
observation [6|, in which the high amplitude oxygens are associated with dark streaks in the 
i?-map, while the low amplitude oxygens are associated with bright spots. Note that the C4 
splitting is characterized by nanoscale domains |6|]. 



VII. LAGRANGIAN FORMALISM 



Here we wish to derive the pairing interaction in the FBM. For this purpose a simple 
approach, akin to the familiar RPA, is to utilise the 1/N expansion technique, where N 
is degeneracy, e.g. vibrational degeneracy n or spin degeneracy n s . In this approach the 
exact solution to the n = 1 quartic vibrator used in the main Ms. will be replaced by the 

n 

(very similar) mean field solution to the n > 1 quartic vibrator [3J. The 1/N expansion 
works well e.g. for the Kondo problem [9], the results remaining physical down to spin 
degeneracy n s = 2. Here we symetrically co-expand in the inverse of the mode degeneracy n 



and the spin degeneracy n s [10|], meaning by expressions such as "1/N" the joint orders 1/n 
and l/n s . The 1/N expansion technique is usually implemented within a Lagrangian/path 
integral formulation lOj, the approach we shall adopt here. We take h = 1, and omit the 
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sum-normalizing factors 1/N%, which can always be restored by inspection. 

Within the usual imaginary time (Euclidean) Lagrangian formulation at finite tempera- 
ture the partition function can be written 

Z = [ VxVce * £(r)dr 



J T>xVce-f° c(T)dT , (25) 



where C (r) is the Lagrangian as a function of imaginary time r, /3 is 1/T (T is temperature, 
taking Boltzmann's constant fee = 1), / T>xT>c implies a path integral over the vibrator and 
fermionic (Grassman) variables. 

The Lagrangian, readily derived from the Hamiltonian (12"U|) . 

£ = £«+4 1 )+4 2) + 4 2 « ) , (26) 
comprises the previously described terms ( JIBl) . ( |T4l) . and ( TTT1) 

4 1} = E - \ E * (* - j) ^ • ( 2? ) 

i,cr 

^ = o E [ mi;2 +a/2 + Xoa; 2 +s/2] 



X0^i+S/2j ) 



2 



5n 

i.a 



,/nn, z — ' 



y 

""i+a/2 7i i+S/2 



A kinetic energy for the fermion variables is introduced on line 1, otherwise the terms are 
the same as in Eq's (1151) . f)14p . and (1171) . In the following we extend the earlier FBM results 
to include the more complete form of coupling derived in Sec. II. 



VIII. RESULTS 

(2) (2) 

We implement a Stratonovic decoupling on the bilinear terms CI and Civ in the La- 
grangian, generating path integrals over new fields z V) i+a/2 and 2 e i+ s/ 2 defined on each 
Cu-O-Cu bond. Later this will make it possible to implement the path integrals over the 
Boson and Fermion fields under the z-path integral. 

Now the partition function is written in terms of the action S 

Z = J VzVxVce- ps , (28) 
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which is again decomposed into the terms 

S = + SP + Sf + S zv + S ze , (29) 



defined by 



S v^ = ? ( U m + W «) X " k -m,a^k,m,a, (30) 



2 

k,m,a 

= (e k — zV TO ) c^ m o .c kimj(T , 

k,m,cr 

C(2) = ^T 2 e e _ 2 V^Ts y> e v 

IK J 1'~ n ' a 1> n ' a y / j — q,— n,a q,n,a; 

q,n,« q,«,n 

5*21) ^ ^ ^— q t — n^^— k,— m,a2-k+q,m+n,ct "I - ^ ^ ^0,0, a ("^ ) ) 

q,k,n,m,a a 

&ze — V ; k,k+q 2: -q,-n,a C k,m,cr C k+q,m+n,o- _ 2 ^o,0,a ( n ) ■ 

q,n,tx,a a 



Here 



x 



(t) = V" e~ iu)nt x n , u n = 2nnT; Bose and z fields, (31) 

n 

c{t) = e~ iUnt c n , v n = (2n + 1) vrT; Fermion fields, 

n 

define the Fourier series converting to Bosonic and Fermionic Matsubara frequencies uj n and 
v ni and a coupling amplitude V'kk' different from that in jsj appears 

i>Z y =4sin[fc Q /2]sin[^/2]. (32) 

The band structure ek and vibrator frequencies ZJ a now involve mean field quantities, an 
effective nearest-neighbor hopping integral t a ( 1331) defined by 

l « = t+z "='-2^M- (33) 

physically meaning that the effective hopping is weakened by the vibrational displacement 
of the oxygen, and vibration frequencies UJ a given by (l34"j) 

mull = X0 - 2z« (34) 

2U / a\ W I 2\ 

= x "-^ n) " + ^ x >•■ 
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physically meaning that the quasiharmonic vibrator is softened by the number of antibonding 
electrons (n a ) a , and stiffened by the effect of the quartic potential. 

On implementing the path integrals over the Boson and Fermion fields under the z-path 
integral (here we drop the a-dependence of the mean field quantities for simplicity), we 
obtain to order (1/iV) the quasi-harmonic action 



& r, [ Z -<l,-n,xi Z -q,-n,y Z -<i,-n,xi Z -q-n,y] A(q, 



n 



q,n 



Z. 



where the matrix A is given by 



A(q, n) 



-nD 2 (n) 
-nD 2 (n) 



2y/nn a 
v 





v 

<\,n,x 

v 

z e 

L *q,n,y 





2,/rm s 



2^/nn s 
v 



f (1 - KR xx (q, n)) -n s R xy (q, n) 

-n s R yx (q,n) % (1 - KR yy (q, n)) 



and the Response functions (RF) are defined by 



D 2 (q, n) 



m?uj 



( uj \ 



coth — + 



<^n,0 



2+4cJ 2 ) V 2 ^; 8TcJsinh 2 (^) 



and (/(e) is the Fermi function) 



D / \ /( e k) — /( e k+q) , a ,3 

JMq, ») = - E £k - £k+q + ^ ^ + X +q ,- 



(35) 



(36) 



(37) 



(38) 



IX. DISCUSSION 



The important regions of A;-space in the Brillouin Zone are the high-density of states 
regions around the saddle points (SP) at k = (0, n) (termed Y), and (n, 0) (termed X). The 
key issue is to understand how the FBM behaves at the SP. Because of the ^-factors, the 
RF are close to being diagonal if we assume only regions around the SP are important in the 
fc-sum. Therefore the A-matrix separates into a component having a pole in the x-channel, 
and a component having a pole in the y-channel. At zero frequency and wavevector, these 
poles correspond to the onset of the phase boundary for the PG at T*. As regards pairing, 
we consider the interaction for a pair to scatter from (k, — k) to (k', — k'), which involves 
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the factor ^ k ,j . This will either scatter where k and k' lie in the neighborhood of the 

point X which involves the factor ^ k ,j ~ 16, in which case the inverse of the A-matrix 
in the ^-channel is involved, or in the neighborhood of the point Y, which involves the 
factor ^0kk') ~ in which case the inverse of the A-matrix in the y-channel is involved. 
The interaction has the classic form required for (i-wave pairing, in that the interaction is 
attractive (in fact close to singular) for (k — k') small, and at large (k — k') the attraction 
dies off (in fact, because of the short-range repulsive interaction U, which is flat in fc-space, 
the interaction actually becomes repulsive at large (k — k')). 

The situation for the full FBM Hamiltonian (|20l) is similar but different in detail for the 
interaction which considers only the X-term in the electronic coupling (fl~9]) {^J , when instead 
of the coupling amplitude ip£ k , we have the amplitude 

X a k y = -2cos(^±^y (39) 

This function leads to a different A- matrix 3| , in which only a component with (i-symmetry 
(the z's in the x- and y- channels having opposite sign) is singular, the other, s-channel, 
being nonsingular. Now scattering around either SP X or SP Y is possible with the the 
factor k , j ~ 4 for either a = x or a = y. However only one combination of channels is 
singular, so we assume that this can be considered as if there is only one effective channel for 
scattering around each SP, just as in the case where both terms in the electronic coupling 
( TT9|) are present. Hence we conclude on the basis of dominance by the SP, that the pairing 
interaction in the case of the full electronic coupling f|T9|) is 4x as large as that in the 
reduced interaction considered earlier. This is the result cited in the manuscript, which 
however needs verification by numerical solution of the gap equation. 
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